function N = reference_basis_2d(basis_type, derivative_type)
switch basis_type
    case "P1"
        N = cell(3,1);
        if (nargin == 1)||(derivative_type == "x")
            N{1} = @(x,y) 1 - x - y;
            N{2} = @(x,y) x;
            N{3} = @(x,y) y;
        elseif derivative_type == "dx"
            N{1} = @(x,y) -1;
            N{2} = @(x,y) 1;
            N{3} = @(x,y) 0;
        elseif derivative_type == "dy"
            N{1} = @(x,y) -1;
            N{2} = @(x,y) 0;
            N{3} = @(x,y) 1;
        end
    case "P1b"
        N = cell(4,1);
        if (nargin == 1)||(derivative_type == "x")
            N{1} = @(x,y) 1 - x - y;
            N{2} = @(x,y) x;
            N{3} = @(x,y) y;
            N{4} = @(x,y) -27*x.*y.*(x + y - 1);
        elseif derivative_type == "dx"
            N{1} = @(x,y) -1;
            N{2} = @(x,y) 1;
            N{3} = @(x,y) 0;
            N{4} = @(x,y) -27*y.*(2*x + y - 1);
        elseif derivative_type == "dy"
            N{1} = @(x,y) -1;
            N{2} = @(x,y) 0;
            N{3} = @(x,y) 1;
            N{4} = @(x,y) -27*x.*(x + 2*y - 1);
        end
    case "P2"
        N = cell(6,1);
        if (nargin == 1)||(derivative_type == "x")
            N{1} = @(x,y) 2*x.^2 + 4*x.*y - 3*x + 2*y.^2 - 3*y + 1;
            N{4} = @(x,y) -4*x.*(x + y - 1);
            N{2} = @(x,y) x.*(2*x - 1);
            N{5} = @(x,y) 4*x.*y;
            N{3} = @(x,y) y.*(2*y - 1);
            N{6} = @(x,y) -4*y.*(x + y - 1);
        elseif derivative_type == "dx"
            N{1} = @(x,y) 4*x + 4*y - 3;
            N{4} = @(x,y) 4 - 4*y - 8*x;
            N{2} = @(x,y) 4*x - 1;
            N{5} = @(x,y) 4*y;
            N{3} = @(x,y) 0;
            N{6} = @(x,y) -4*y;
        elseif derivative_type == "dy"
            N{1} = @(x,y) 4*x + 4*y - 3;
            N{4} = @(x,y) -4*x;
            N{2} = @(x,y) 0;
            N{5} = @(x,y) 4*x;
            N{3} = @(x,y) 4*y - 1;
            N{6} = @(x,y) 4 - 8*y - 4*x;
        end
    otherwise
        error("Invalid basis type.");
end
end